Statistical properties of the linear tidal shear 
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Large-scale structures originate from coherent motions induced by inhomogeneities in the primeval 
gravitational potential. Here, we investigate the two-point statistics of the second derivative of the 
potential, the tidal shear, under the assumption of Gaussianity. We derive an exact closed form 
expression for the angular averaged, two-point distribution of the shear components which is valid 
for an arbitrary Lagrangian separation. This result is used to write down the two-point statistics 
of the shear eigenvalues in compact form. Next, we examine the large-scale asymptotics of the 
correlation of the shear eigenvalues and the alignment of the principal axes. The analytic results 
are in good agreement with measurements obtained from random realizations of the gravitational 
potential. Finally, we show that a number of two-point distributions of the shear eigenvalues are 
well approximated by Gaussian bivariates over a wide range of separation and smoothing scales. 
We speculate that the Gaussian approximation also holds for multiple point distributions of the 
shear eigenvalues. It is hoped that these results will be relevant for studies aimed at describing the 
properties of the (evolved) matter distribution in terms of the statistics of the primordial shear field. 
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I. INTRODUCTION 



In the currently favored ACDM cosmology, the galax- 
ies and the large-scale structures we observe today are 
thought to arise from the hierarchical growth of initially 
tiny Gaussian fluctuations. Galaxy surveys [ij reveal that 
large-scale structures on scale '--^ 10— 100 /i^^Mpc form a 
predominantly filamentary network whose principal con- 
stituents - clusters, filaments and walls - accumulate on 
the boundary of large voids. Numerical simulations and 
semianalytic approaches have been very successful in re- 
producing the observed matter distribution [3], whereas 
analytic models based on the spherical infall [31 predict 
mass functions, merging histories and spatial clustering 
of bound objects that are in reasonable agreement with 
the observations [1]. In parallel, several methods have 
been proposed to describe quantitatively the structures 
observed both in the Universe and in the simulations. 
The full hierarchy of correlation functions, the two-point 
statistics in particular remains the most widely used 
statistical tools to distinguish between different scenar- 
ios of structure formation and constrain the cosmological 
parameters. Topological estimators such as Minkowski 
functionals [6*1 provide useful complementary informa- 
tion on the morphological characteristics of the filamen- 
tary network. Furthermore, various identification algo- 
rithms abstracting the spatial patterns in points, lines, 
etc. have been proposed in an attempt to improve upon 
current topological measures [7]. 

Although the sequence in which large-scale structures 
form is still a matter of debate, many lines of evidence 
suggest that the filamentary pattern seen in observations 
and in N-body simulations is a consequence of the spatial 
coherence of the initial tidal shear ^] . While the spheri- 
cal infall model [3] captures the essential features of grav- 
itationally induced collapse, the primeval shear field has 



also been shown to play a crucial role in the formation of 
nonlinear structures 0] . As demonstrated in [l3| , the in- 
clusion of nonsphericity in the collapse dynamics yields a 
better fit to the halo mass functions measured in N-body 
simulations. Yet another important manifestation of the 
tidal shear is the alignment of shape and angular mo- 
mentum of objects [11, 12, 13, 14, 15] . Numerical stud- 
ies of the ACDM cosmology report strong correlations 
in the alignment of g alaxies, haloes, massive clusters, or 
voids [16, ITI. [TsI. [191 [20| . reflecting the coherence of the 
matter distribution out to large distances. 

In the "Cosmic Web" picture outlined in [1], the cor- 
respondence between structures in the evolved density 
field and local properties of the linear tidal shear should, 
in principle, allow us to estimate the morphology of the 
matter distribution. In practice however, this corre- 
spondence has not been much exploited principally be- 
cause of the lack of theoretical results. In spite of the 
progres s made in the analysis of Gaussian random fields 
|2ll |22|. [23j applied to the formation of large-scale struc- 
tures, the statistics of the shear has received little at- 



tention. Doroshkevich [2j| first calculated the probabil- 
ity distribution of the shear eigenvalues and ascertained 
the amount of material being incorporated in sheetlike 
structures or pancakes. Reference [25| reexamined the 
formation of these pancakes and derived a distribution 
function for the largest eigenvalue of the shear tensor. 
Reference '^Sj computed conditional probability distri- 
butions for individual shear eigenvalues and obtained an 
analytic approximation to the halo mass function. Also, 
[131 explored the two-point correlation of the tidal shear 
components, but they did not discuss probability distri- 
butions. 

In this paper, we carry out the analysis of the 2- 
point statistics of the linear tidal shear at two distinct 
(Lagrangian) positions. We extend the study of [28j . 
who derived an expression for the shear 2-point statis- 
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tics smoothed at two different scales, but evaluated at 
a single position. This work is essentially intended to 
provide theoretical results that could improve the statis- 
tical description of the tidal shear and, therefore, of the 
Cosmic Web. Section Hill is devoted to the calculation of 
the 2-point distributions. A careful examination of the 
small and large scale behavior of the 2-point distribution 
of shear components suggests a compact expression for 
the joint distribution of shear eigenvalues. This result is 
used in Sec. lIVI to explore the asymptotic behavior of the 
correlations of shear eigenvalues and principal axes. In 
Sec. |Vl conditional 2-point distributions obtained from 
random realization of the linear shear field are compared 
with theoretical Gaussian distributions. It is argued that 
Gaussian multivariates provide a good description of the 
n-point distributions of shear eigenvalues at all separa- 
tion and smoothing length. This suggests the possibility 
of implementing nonspherical collapse in current analytic 
models of structure formation using well-known results of 
the theory of Gaussian random fields. 



II. SHEAR 

The comoving Eulerian position of a particle can be 
generally expressed as a mapping 



x= q+S{q,a) , 



(1) 



where q is the Lagrangian (initial) position, a is the scale 
factor and S is the displacement field. Continuity im- 
plies that the Eulerian density contrast d{x,a) is given 
by the reciprocal of the Jacobian for the transformation 
([H), l-|-^(a;, a) = \det (dq/dx) \. Singularities occur in 
this mapping whenever at least one of the eigenvalues is 
positive, signaling crossing of particle trajectories at that 
Eulerian point. The initial deformation tensor or strain 
field Dij = diSj {di = d/dqi) thus encodes important 
information on the collapse of large-scale structures. 

In the Zeldovich approximation, the displacement 
field is S{q,a) = -D+(a)V$(g), where $(g) = 
(l){q, a) /4TrGp,nio-)a'^ D^{a) is the perturbation potential, 
(f){q, a) is the Newtonian gravitational potential, pm is the 
average matter density and D-\-{a) is the linear growth 
factor (normalized such that D+{a) = 1) [5|. The strain 
tensor now is proportional to the second derivatives of 
the perturbation potential. For convenience, we will in- 
troduce the real, symmetric tensor 



cr dqidqj 



(2) 



where a is the root-mean-square (rms) variance of den- 
sity fluctuations 6{q) = A$(q) linearly extrapolated to 
the present time. We shall henceforth refer to as the 
shear tensor. Notice that is dimensionless while $(q) 
has the unit of (length) 2. We will also assume that these 
fields are smoothed at some characteristic scale Rs with 
a spherically symmetric window W{r, Rs). Although 



many choices are possible for such a filtering function, 
we will adopt a top-hat filter throughout this paper, so 
that the variances are related to spherical volumes of ra- 
dius Rs- 

Let A = diag(Ai, A2, A3) be the diagonalized shear ten- 
sor. For Gaussian initial conditions, the 1-point proba- 
bility distribution of the shear eigenvalues derived in [23| 
can be written as 



P(Ai,A2,A3) 



15^ 



g-Qi(A,A)^(^) , (3) 



87rV5 

where, for shorthand convenience, 



Qi (X, Y) = - [5tr (XY) - (trX) (trY)] (4) 

is some (indefinite) quadratic form over the space of real 
matrices, and 



A(A) = dct {X^-')=1[{K-Xj) 



(5) 



is the Vandermonde determinant in the arguments 
•^11-^2, A3. Our definition ([2]) of the shear tensor makes 
the probability P(Ai, A2, A3) independent of the filtering 
scale Rs- For instance, one finds that, for ambient field 
points, the probability of all three eigenvalues being pos- 
itive is P{+ + +) = 0.08, and that of the configurations 

(+ H — ) and (-1 ) is 0.84. Note, however, that these 

values depend strongly on the peak height, v = S/a, of 
the region under consideration, the highest density peaks 
being predominantly spherical [2^ l30j|. 



III. TWO-POINT STATISTICS 

Desjacques [2^ extended the results of [l^l to the joint 
statistics of the shear smoothed at different scales. How- 
ever, he confined his calculation to the case in which the 
joint distributions are evaluated at a single Lagrangian 
position. Here we address the general case and calculate 
the joint distribution of the shear components (9i) ^^^id 
iki{(l2) for arbitrary Lagrangian separations r= 52 — 9i- 
We shall assume throughout this paper that the ini- 
tial fluctuations are described by the Gaussian statistics. 
This assumption is remarkably well supported by the lat- 
est measurements of the cosmic microwave background 

(CMB) [m. 



A. Shear correlations 

We take the components Cij(9i) and £,ki{Q2) to be 
smoothed at two different (comoving) scales Ri and R2, 
respectively. The spectral parameter 
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a'iO'2 Jo 



dlnfc A^(fc) W{k, Ri)W{k, R2) , (6) 
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< 7 ^ Ij is a measure of the correlation between these 
scales. Here, Aj{k) = PPs{k)/2TT'^ is the dimension- 
less power spectrum of the density field, W{k, Ri) is the 
Fourier transform of W{r^ Ri), and ai is the rms variance 
of density fluctuations 6{q) smoothed at scale Ri. 

Statistical isotropy and symmetry imply that, in posi- 
tion space, the 2-point correlation functions of an arbi- 
trary symmetric tensor field S,ij{q) must be of the form 

ij ^Ira 

) , (7) 

where r = 1 5i — | , h — ri / r is a unit vector and the 
functions 'l'i(r) depend on r only. This is the most gen- 
eral ansatz for the isotropic sector of the fourth order 
correlation function {^ij^im){r). Symmetry requires that 
f appears in even number pairs. In the case of a scalar 
(spin-0) tensor such as the linear tidal shear defined in 
Eq. (121), *2 = *3 and *4 = ^s- Notice that Eq. ^ 
holds regardless of the statistical properties of the gravi- 
tational potential. However, when is Gaussian, the 
functions ^'^ may be conveniently expressed as 



dlnfc A2(fc)j4(fcr 
dlnfcA2(fc) 
dlnfcA2(fc) 



(8) 
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where A^(fc) = A^(fc)W^ilV2/(o'icr2) and je{x) are spher- 
ical Bessel functions of the first kind. The can be 
equivalently expressed in terms of the auxiliary functions 



Jn = nr " /Jds 



,n-l 



27|,|3a,l33(, where 



^(r) = *i(r) + 10*3(r) + IS^sW 



(9) 



is the cross correlation between the density enhancement 
6/(7 smoothed at two different scales. In the limit r ^ 0, 
both ^E"! and ^"3 vanish while ^5 tends towards 7/15, so 
that ■0^7. 



B. Two-point probability distribution 

Owing to the symmetry of ^ij , only six components of 
the shear are independent. Following the notation of [1^, 
let ^ = {£_A, ^ = 1, . . . , 6} designate the six-dimensional 
vector whose components are equal to the components 
ij = 11, 22, 33, 12, 13, 23 of the shear tensor. The joint 
probability distribution P (^1,1^2;'") of the shear tensor 

= ^ijili) and ^2 — ^ij{(l2) is given by a multivariate 
Gaussian whose covariance matrix C has 12 dimensions. 
This 12 X 12 matrix may be partitioned into four 6 x 
6 block matrices, Ci = in the top left corner. 



C2 = (^2^J) in the bottom right corner, B = (^i^J) and 
its transpose in the bottom left and top right corners, 
respectively, where 

(10) 

and I is the 3x3 identity matrix. 

Unlike Ci and C2, the cross correlation matrix B gener- 
ally is a function of the separation r. Using the harmonic 
decomposition of the tensor products f (g) • • • ® f which 
appear in Eq. ([7]), B(r) can be written as follows : 



B(r) = 



Bi(r) B3(r) 



B3(r) B2(r) ; ' 
with 3x3 block matrices 

Bi(r) = -L^(.)Mi+ ^ Bt™(r)r;"(r) 



(11) 



i=2A 



B2(r-) = ^V(0I+ E ^'rir)Yr{r) 

e=2,4 

£=2,4 



(12) 



Y™{r) are spherical harmonics and Bj'™(r) are 3x3 ma- 
trices which satisfy (B^''")t — (— 1)™B^'™ (Bi(r) are real- 
valued matrix). An explicit calculation of these matrices 
is not necessary as we will focus on the contribution of 
the monopole terms. Again, symmetry implies that only 
the harmonics with multipoles £ = 0, 2, 4 and even m 
appear in the decomposition. Furthermore, it is worth 
noticing that the joint probability density P(^i,^2; r) is 
invariant under any arbitrary rotation of the coordinate 
system [s^ . 

Piei,e2;r') = P(RaR^,R6R^,R^r') - P{^i,^2,r) , 

(13) 

where primes denote quantities in the new coordinate 
frame. However, in a given coordinate system, its func- 
tional form changes when r moves over the unit sphere. 
Namely, the transverse and longitudinal components of 
the 2-point shear correlation vary with the orientation of 
the separation vector. Therefore, the anisotropic struc- 
ture of {^ij^im)ir) should come as no surprise [27| . 



C. Angle average distribution 

We are primarily interested in the angular average of 
the 2-point probability distribution of the linear tidal 
shear, 



P(ei,6;0 = ^ Jdn,P{^ui2;r) 



(14) 



which is a function of the separation r solely. To get in- 
sight on the functional form of P(^i, S.2]f)i we will exam- 
ine the small-scale {r <^ 1) and large-distance asymptotic 
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(r 1) behavior of the probability density P(^i,^2,'")- 
In both regimes, the cross correlation matrix B can be 
expressed as B(r) = B + dB{r), where 5B is a small per- 
turbation and B is the zero order contribution, which is 
either B — 7C1 (when r — > 0) or B = (when r —^ 00). 
The quadratic form which appears in the probability dis- 
tribution P(fi,C2; r), 



1 



i2ny |detC|i/2 



.-Q(«l,«2;r) 



(15) 



can be computed easily using Schur's identities. Here, 
detC is the determinant of the covariance matrix C. Ex- 
panding the exponential in the perturbation 5B and av- 
eraging over directions f, we find 



^ ydr!,exp[-Q(6,6;r')] 

(■0 - 7) 



1 - 2- 



1-72 



1-7^ 



(Qi(a,Ci) + Qi(6,6)) 



-[Ql(«l:?l)+Ql(?2,?2)-27Ql(?l,«2)]/(l-7') 



-^1(^1,6)]} 

(r « 1) , (16) 

« [l + 20Qi(ei,6)] e-«i(«-«i'-«i(«-«^) 

(r » 1) , (17) 

to first order in SB. Interestingly, these perturbative ex- 
pansions precisely match the small-scale and asymptotic 
large r behavior of the function exp[— (52(^1 : r)], where 



Qi(a,ei) + Qi(6,6)-2V'gi(a,6) 
1 -1/.2 



^2(^1, 6; 

1 — 'UJ" 

(18) 

This strongly suggests that the 2-point probability dis- 
tribution P{^i,£.2','r) of the shear tensor may be written 
exactly as 



1 
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20 V 27r 



(19) 



This distribution depends on the separation r through 
the density correlation ip{r) only. Gaussianity and in- 
variance under rotation requires that P be a function 
of ti{^f), (tr^,)2, tr(CiC2) and tr^itr^2 solely. Although 
we have not been able to rigorously prove that Eq. ([TO]) 
correctly describes the 2-point distribution in the inter- 
mediate region, we have found that it agrees with the 
result of a direct numerical integration of Eq. (fTil) (for 
various choices of ^1, ^2 and r) up to the numerical ac- 
curacy. Furthermore, the fact that, in the limit r <g; 1, 
Eg . (1191) reduces to the joint probability density derived 
in 1281 is another indication of correctness. 



D. Joint distribution of the eigenvalues 

We now choose a coordinate system such that the 
coordinate axes are aligned with the principal axes of 
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FIG. 1: The density correlation 'ip{r) and the parameter /3(r) 
(cf. text) as a function of the ratio r/Rs for two different 
smoothing lengths Rs = 0.1 and 1 /i~^Mpc. 



^1. Let Ai and A2 be the diagonal matrices consist- 
ing of the three ordered eigenvalues xi > X2 > X3 and 
yi > 2/2 > 2/3 of the deformation tensors ^1 and ^2, re- 
spectively. The principal axes are now labeled according 
to this ordering. With this choice of coordinate, ^1 = Ai 
and ^2 = RA2R^, where R is an orthogonal matrix that 
defines the orientation of the eigenvectors of ^2 relative 
to those of ^1. To preserve the orientation of the prin- 
cipal axis frames, we further impose the condition that 
the determinant of R must be -1-1. Namely, R belongs 
to the special orthogonal group SO (3). The properties 
of the trace imply that Qi{£,2,^2) = Qi{^2,^2), while 
the term Qi{£,i,^2) = Qi(Ai,RA2R^) depends on the 
rotation matrix. 

To obtain the joint probability distribution of the or- 
dered eigenvalues of the shear tensor, "angular" vari- 
ables, such as those appearing in Qi{S.i, ^2), have to be 
integrated over. The volume measure d^ for the space of 
real 3x3 symmetric matrices can be expressed in terms 
of the nonincreasing sequence of eigenvalues Ai (= Xi or 
Vi) as 



d^ = 87r2 A(A)d3AdR 



(20) 



Here, dR is the Haar measure on the group SO (3) nor- 
malized to / dR = 1, d'^A = dAidA2dA3 and A(A) is the 
Vandermonde determinant equation ([5]) . When the rota- 
tion matrices R are parametrized in terms of the Euler 
angles < tpjip < 27r, < i9 < tt, the Haar measure takes 
the familiar form 



dR ~ sin d dipdddip 
87r-^ 



(21) 



Since the quadratic form Q depends only on the relative 
orientation of the eigenvector triads of ^1 and ^2, we can 
immediately integrate over one of the SO (3) manifolds. 
The relevant volume is Stt'^ /A = 2Tr^ . The factor 4 comes 
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from not caring whether the rotated axis points in the 
positive or negative direction [2^ . The essential problem 
is the calculation of the integral over the manifold that 
defines relative, distinct triad orientations, 

dRexp[/3tr(RA2R^Ai)] , (22) 

/SO(3) 

where we have defined — (15/2)'(/'/(l — ■0^)- 

There is no analytic, closed-form solution to this in- 
tegral, although an exact determinantal expression was 
derived when averaging over the unitary group U(N) [H, 
IsgI I . An asymptotic representation can be obtained when 
/3 » 1. For reasonable values of Rs, this occurs when 
the separation r is less than a few smoothing radii (see 
Fig. [1]) . In general, the integral ((22)) can be expressed as 
a hypergeometric series with the symmetric 3x3 matri- 
ces /3Ai and A2 as argument (see Appendix In the 

notation of [IJ, ^ = o-fo^'^'' where, at second order in (3, 
oF^^\pA,,A2) = 1 + ^trAi trAa + ^(trAi)2(trA2)^ 



90 



[3tr(A2) - (trAi)2] [3tr(A2) - (trAa)^] .(23) 



Higher order terms are intentionally omitted as they are 
not used in the present analysis. 

A straightforward calculation shows that the eigen- 
value joint probability distribution P(x, y; r) evaluates 
to 



(24) 



where 



Q12 = 



Qi(Ai, Ai) + Qi(A2, A2) + IV^(trAi) (trA2) 
(1-^2) 



(25) 

Notice that, in the limit ^ 1, the joint probability 
distribution P{x,y;r) tends, as it should be, toward the 
product of the individual 1-point probability distribution. 
Using Bayes' theorem, we easily derive a conditional dis- 
tribution P{x\y;r) for the shear eigenvalues Xi given the 
UiS and a separation r, 

P{x\y;r) = ^ (l - V^)-^ ^i^f (/JAi 



where the quadratic form Qi\2{Ai, A2;r) is 



A2)e- 



-Qi 



A{x) 
(26) 



gi(Ai, Ai) + gi(A2, A2) + |V;(trAi)(trA2) 



Ql\2 

± — 

(27) 

A direct numerical integration convinced us that the 
probability distribution (|26p is normalized to unity (we 
used the multidimensional integrator DCUHRE described 
in [3^). We believe that Eqs. ^ and ^ are exact 
expressions for the 2-point and conditional probability 
distribution function of the shear eigenvalues. They gen- 
eralize the results obtained in 123, 1231 . 



IV. ASYMPTOTICS OF CORRELATION 
FUNCTIONS 

Instead of attempting a brute force calculation of the 
correlation functions through a direct integration of the 
probability density (P^ . we will examine the large-scale 
asymptotic behavior solely. We will nonetheless infer an- 
alytic approximations to the correlations of shear eigen- 
values which are accurate on all scales. 



A. Eigenvalues 

In order to derive the correlation function for the eigen- 
values of the shear tensor in the asymptotic regime r ^ 1 
for which ip <^ 1 and (3{ip) ^ 1, we transform to the new 
set of variables {vi, ei,Pi, i — 1, 2}, where vi = xi + X2 + 
X3, ei — (xi ~ x:i)/2i>i and pi = (xi — 2x2 + X'i)/2vi. 
The variables (j^2, 62,^2) are defined as a function of the 
eigenvalues Ui in a similar way. and Pi are the shear 
ellipticity and prolateness, respectively. Our choice of or- 
dering impose the constraints > and — < Pi < e^. 

The cross correlation function C,ij{r), or 2-point con- 
nected moment of the shear eigenvalues Xi and yj is de- 
fined as 



{xi){yj) +C,ij{r) ^ d xd yP{x,y;r) Xiyj 



(28) 



The integration over the variables pi and is straight- 
forward to second order in (3. In this respect, notice 
that the volume measure d'^x and the Vandermonde 
determinant A{x) are d^a; = (2/3)j^i di^ideidpi and 
A(a;) = 2iyf ei (e^ — Pi), respectively. Furthermore, with 



the help of the series expansion (|23|) about (3 = 0, the 
average over the relative orientation can be reduced to 



0^0 



(3) 



P P2 2 

'gi^lI^2 + Yl^ll'2 



^-{3el+pl) {3el+pl) 
(29)^ 

in the aforementioned coordinates. The integration over 
the peak heights vi and 1^2 is then easily performed and 
gives 

Ur) = l^ir) + ^{x.ri^'ir) (30) 

at second order for all three eigenvalues. To derive this 
result, we have used the following expressions for the 
mean of the shear eigenvalues [U, [2^ , 



{xi) = {xs) = 



{X2) = 



(31) 



Hence, there is no second order contribution to C22if)- 
It is worth remembering that the variance (xf) of each 
eigenvalue is 



(4) 



137r - 27 
307r ' 
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15 



(32) 
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FIG. 2: A comparison between a measurement of the auto- 
and cross- correlation of the shear eigenvalues from realiza- 
tions of Gaussian random fields (symbols) and our analytic 
approximations (curves). All the correlations but ^33 have 
been shifted vertically for clarity. On large scale, the corre- 
lations asymptote to the exact result i/)(r)/9 shown as the 
dotted curve. On small scale, however, Cjj('^) defined 
in Egs. 1301 and I36j provides a better fit to the measurement 
(dashed curve) . The bottom panel shows the fractional error. 
Notice that the density correlation ip{r) used for the compar- 
ison is calculated from the numerical realizations to account 
for the missing power at small and large scales. 



close to the value of 1/9. This suggests that the corre- 
lation functions Cii('") £^re well approximated by tp{r)/9 
on all scales. This is not entirely surprising since, in 
the case where the correlations between eigenvalues are 
equal, the constraint J^ijdji''') — V'('') would imply 
C,,(r)=V(r)/9. 

In order to verify this assumption, we have gener- 
ated random realizations of the potential field ^{q) for 
the ACDM cosmology considered here on a 256^ mesh 
of size 250 /i~^Mpc. The eigenvalues of the shear ten- 
sor are computed on the same grid using standard FFT 
(Fast Fourier Transform) techniques. More precisely, 
the Fourier modes of the shear are computed using the 
relation ^ij{k) = kikj^{k)/a. Once (,ij{k) is Fourier- 
transformed back, the shear eigenvalues well as the 
density v — xi + X2 + X3 are calculated at each grid 
point. Lastly, after having checked that the measured 
variances (xf) match well the analytic expectation ([5^ . 
we calculate correlation functions of the shear eigenval- 
ues. Note that the gravitational potential is smoothed 
on scale Rs = 1 h~^Mpc. 

In Fig. [21 the correlations (22 {r) and C33 {r) (recall that 
Cii = C33) averaged over the realizations are shown as 
empty symbols. These measurements are compared to 



the asymptotic scaling (pO)) and to the following analytic 
estimates 



CiiM 

C22(r) 



i-0(r) -I- —ib^M , 



(33) 



which are designed to asymptote to the variances given 
in Eq. ([5^ . The density correlation function ip mea- 
sured from the simulations is used for the evaluation of 
those theoretical expressions. Figure [2] clearly demon- 
strates that, while the linear approximation -0/9 is in ex- 
cellent agreement with the measurements in the asymp- 
totic regime, it deviates at least 10 percent at small dis- 
tance, r <3 /i^^Mpc. By contrast, Cm('') as defined above 
achieves a fractional error no larger than 2 percent for 
separations less than ~ 30 h~^Mpc. 

These results are readily extended to the cross corre- 
lations Cij, i ^ j- Proceeding in a similar way, Cii(^) can 
be rearranged as 

Q,{r)=^-^{r) + ^{x,){y,)i,^{r) (34) 

to second order in tp{r). Estimating the cross correlations 
at zero lag Cij(O) = i^i^j) from the theoretical proba- 
bility distributions P{xi,Xj) proves difficult. A numeri- 
cal integration gives the following hypothesized rational 
forms 



{X1X2) 



{X2X3) = ^, 



{X1X3) 



27 - 677 
SOtt 



(35) 



for which the constraint {{xi +X2 -l-xa)^) = 1 is satisfied. 
This motivates the interpolation formulae 



1 



Ci2(r) = C23(r) = - Hr) ~ - ^\r) , (36) 



which match reasonably well the large- and small-scale 
behavior of (Fig. [2|) . Finite grid size effects may be 
responsible for the slight offset (roughly 2 percent) of the 
cross correlation (^13 relative to the theoretical prediction. 

Numerical investigations indicate that dark matter 
haloes do not form randomly in the initial conditions, 
but rather preferentially close to the peaks of the density 
field [13] ■ To assess the extent to which biasing affects 
our results, we adopt in a first approximation the usual 
critical density criterion issued from the spherical infall 
model. As first recognized in [4l|, the correlation func- 
tion of regions lying above a certain density threshold v is 
enhanced relative to that of the density correlation ilj(r). 
Likewise, the correlations of shear eigenvalues restricted 
to regions with density larger than a given threshold v 
are also amplified. On large scales and in the regime 
I' 3> 1, we find 



C,^^{r\ > V) 



1 /IN 
- + v{Xi\v) 



■)p{r) , 



(37) 
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conventional to use 



12 3 4 

peak height v 

FIG. 3: Correction factors 6ii(A, v) (see text) as a function of 
the threshold height u when the alignment of principal axes 
is restricted to those regions where all shear eigenvalues are 
positive. The bottom panel shows the difference &ii — 633, 
which demonstrates that the alignment between major axes 
is slightly stronger than the correlation of minor axes. The 
correlation of the intermediate axis is strongly suppressed over 
the whole range of peak height. 



where the conditional average eigenvalue 



3 ' 



(38) 



depends linearly on the peak height i^. This should be 
compared to the correlation function of thresholded re- 
gions, which is ^(r| > i/) w 1/ "*■)/) (r) in the same limit [4]| . 
Hence, CaiA > ^) ^Iso exhibits the usual linear amplifi- 
cation factor v'^ of dense regions. It would be interesting 
to estimate the extent to which the large-scale bias varies 
when constraints are imposed on the shear eigenvalues. 
This calculation is postponed to a future work. 



B. Principal axes 

We now turn to the correlation of the shear princi- 
pal axes. Since it is computationally expensive to mea- 
sure such a correlation from numerical realizations (di- 
rect summation must be employed), we will only present 
an analytic estimate which is valid on large scales. We 
also discuss its dependence on the peak height. Note 
that similar calculations for the correlation of angular 
momentum Li oc Cijk£,jilik, where Iik is the inertia ten- 
sor of some Lagrangian region, can be found in [s^, 
for instance. 

Let hi designate the unit vector in the direction of the 
major, intermediate, or minor axis of the shear. It is 



(39) 



as a measure of the alignment between principal axes |42 . 
l43j . so that rjijir) = in the absence of any correlation. 

We will only consider the correlations -qn since the 
calculation of r/y proceeds along similar lines. We 
parametrise the rotation matrix R in terms of the Eu- 
ler angles < '0 < 27r and Q < •& < tt. We adopt the 
XYX, YZY, and ZYZ convention when working out the 
correlation of major, intermediate and minor axes, re- 
spectively, so that cos d always is the cosine of the angle 
between the considered axes. The average is performed 
over the independent components of the shear tensor. 



cos^i?-!) P(a,6;r-) , (40) 



where, in the limit r 3> 1, the 2-point probability dis- 
tribution P{^i,£,2', r) reduces to Eq. (|17p . The quadratic 
form (3i(Ai, RA2R^) can be expanded in terms of the 
Wigner D-functions Pj^^ , I being the index of the rep- 
resentation. These 3D harmonics generate irreducible 
representations of the three-dimensional rotation group 
and, therefore, form a complete orthogonal set of func- 
tions defined on SO (3) itself. Invariance under reflec- 
tions implies that the quadrupole rotation matrices with 
TO = 0, ±2 appear in the harmonic decomposition of the 
quadratic form Qi as follows : 

gi(Ai,RA2R^) = itrAitrAa 



15e_ 



^1 25 



+2 ^1^2 



(41) 



The explicit form of these harmonics in the ZYZ represen- 
tation is given in Table IH It is worth emphasizing that, 
apart from the "traceful" contribution l/2trAitrA2, 
Qi(Ai,RA2R^) depends on the three "shape" param- 
eters e_, ei, and 62 solely, because points on SO (3) 
truly have only 3 degrees of freedom. These param- 
eters are given by e_ = (1/3) (a:;i3 + a;23)(?;i3 + 2/23), 
ei = (xi - a;2)/(xi3 a;23), £2 {vi ~ 2/2)/(2/i3 + 2/23), 
where Xij = Xi — Xj and yij = yi — yj. The alignment 
is simply rjij = 2/3 {D'^^) with our parametrization of 
the rotation matrix. The addition of angular momentum 
then yields 



dR 

SO(3) 



COS^ 



gi(Ai,RA2R^) 



l6 



(42) 



after averaging over the angular variables. Lastly, the 
integral over the volume measure d'^xd'^y is easily per- 
formed in the coordinate system (i/, e,p), where 



_ J \i>iV2 (3ei — pi) (3e2 — P2) minor axis 
3i^ii^2 (3ei -f pi) (3e2 -I-P2) major axis 



(43) 
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FIG. 4: A comparison between the conditional 2-point probability distribution P{yi\xi;r) measured from realizations of Gaus- 
sian random fields and the Gaussian approximation equation (I51|) . Results are presented for the largest and intermediate 
eigenvalues only (left and right panels respectively). Symbols show the measurements of the conditional probability for a sep- 
aration 2 < r < 2.5 /i~^Mpc and for Xi = (xi) ± 2ai. The correlation used in the analytic estimate shown as the dashed 
curve is directly calculated from the random realizations of the shear field. The bottom panels show the fractional error. 



Let us choose vi — 1^2 — i' for illustration purposes, and 
perform the integration over the domain defined by > 
and \pi\ < Ci {i = 1,2). We find that, in the large-scale 
limit r 3> 1, the alignment r]ii{r\ > v) of thresholded 
regions evaluates to 

Vuir\>iy)^^^ir)^0A3i^ir) , (44) 

regardless of the peak height and the axis under consider- 
ation (major and minor). This is unsurprising given the 
invariance of the integral over the asymmetry parameters 
under the reflection pi —pi. For the intermediate axis, 
a similar calculation yields e_ = (4/3)z^ii^2PiP2- This im- 
plies 7722 (''I > J^) = at leading order since the 1-point 
probability P(e,p|i/), 

Pie,pW) = }^,^e{e'-p')e-^'^'('^'+^') , (45) 

is symmetric about p — 0. 

Unlike the conditional correlation Cm(^| > 1^) of shear 
eigenvalues, the alignment r]ii{r\ > v) between the shear 
principal axes is insensitive to the peak height. How- 
ever, restricting the domain of integration to the region 
where all shear eigenvalues are positive, for instance, 
can introduce a dependence on the threshold height. 
Such a constraint naturally arises in models of struc- 
ture formation to characterize the Lagrangian regions 
which collapse into dark matter haloes. The domain 
where the lowest eigenvalue is positive corresponds to 
the interior of the triangle bounded by {e,p) = (0,0), 



(i,— i), and (^,5). The conditional correlation can 
be conveniently expressed as bii{A,h')r]{r\ > i'), where 
ri{r\ > v) = 27/207r?/;(r) and &ii(A, v) is a correction fac- 
tor resulting from the restriction to the triangular domain 
A3 > 0. After some manipulation, we find 




These expressions are plotted in the upper panel of Fig. [3] 
as a function of the threshold height. Clearly, the align- 
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TABLE I: Quadrupole Wigner D-functions ©^^^ {(p,'d,tp) in the ZYZ representation. Harmonics with m\,m2 = ±1 are not 
shown since they are not needed for the present analysis. 







1712 = —2 m2 = 


m2 = 2 


mi 


= -2 


i(l + cos^)'e2^*'+2*'^ y^fsin^^e^^*' 


i(l-cosi?)'e2'^-2''^ 


mi 


= 


yi'sin^^e^^'* i(3cos2,9-l) 


sin^^e-^'* 


mi 


= 2 


i (1 - COS^)2 ^~2^^+2^i, ^3 sin2^ ^-2-^ 


i(l + cos e~2'^-2''/' 



ment between major and minor axes is strongly sup- 
pressed when V <2. The correction factors 6ii(A, v) are 
always less than one, to which they asymptote in the 
limit of large threshold v. Furthermore, in the range 
1/ ^ 1 — 3, we have bn >b^'i so that the correlation be- 
tween major axes is stronger by a few per cent. Shown in 
the lower panel of Fig. [3] is the difference 611 — 633. It is 
maximal when the integral of the probability distribution 
P{e,p\v) over the triangle is about one-half. This occurs 
when the mean ellipticity (e|i^) = i / {s/ roughly is 
1/4— 1/2, i.e. when v k, 2. Restricting the integration 
domain to A3 > thus induces a nonzero, albeit small, 
correlation between intermediate axes. ri22{f'\ > v) does 
not exceed ~ O.O1'0, a value reached when v ^ 0.5. We 
emphasize that the alignment between the principal axes 
of the shear field is in all cases proportional to the den- 
sity correlation on large scales. This is consistent with 
the findings of [33, 44]. 

We have shown that the initial alignment of the shear 
principal axes is significant only for those regions which 
collapse into haloes of mass M >Mi,. We will, nev- 
ertheless, not extend the discussion to the statistics of 
dark matter halo and galaxy shapes as it is unclear to 
which extent the correlations detected in N-body sim- 
ulations of CDM cosmologies [l3|, or in galaxy surveys 
(l6i] . reflect the large-scale alignment of the initial shear. 
Nonlinear effects, such as anisotropic accretion or re- 
laxation following collapse, could plausibly enhance or 
erase the large-scale coherence of the primordial tidal 
field [3, [4^. Indeed, even in the artificial case where 
all dark matter haloes are perfectly aligned would the 
alignment and clustering of galaxies be negligibly af- 
fected (45|. We conclude by noticing that primordial 
non-Gaussianities characterized by a local mapping of 
the form 0ng = - /nl(0' - (0')) [ii, IH Ef would 
not affect the correlation of principal axes. 



V. A GAUSSIAN APPROXIMATION 
A. Shear eigenvalues 

Explicit expressions for the 1-point probability distri- 
bution P{xi) of the individual shear eigenvalues Xi can 
be found in [H, [2^ . It is worth noticing that , although 
the variables Xi{q) are not Gaussian random fields, their 
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FIG. 5: The conditional probability density P{yi\xi;0) eval- 
uated at a single Lagrangian position but at two different 
smoothing lengths Ri and 7?2. The triangles, squares, and 
crosses show P{yi\xi;0) measured from random realizations 
of the potential smoothed on scale R2 ~ 2.5, 3, and 5 /i~^Mpc, 
respectively, while keeping _Ri = 2 h'^Mpc fixed. The dot- 
ted and dashed curves indicate the Gaussian approximation 
when xi — (xi) + 2ai and (xi) — 2a\, respectively. The bot- 
tom panel shows the fractional error. 



1-point probability distributions are very close to Gaus- 
sian. More precisely, P{x2) is indeed (fortuitously) ex- 
actly Gaussian, whereas the probability densities P{xi) 
and P{x3) show a small positive skewness. 



skewness — 



33/2 (54 - 1 77r) 
(137r-27) 



3/2 



0.060 , 



(49) 



which refiects the fact that the large and small tail of 
these distributions, respectively, is slightly more pro- 
nounced. 

To assess the extent to which 2-point statistics of 
the shear eigenvalues deviate from Gaussianity, we have 
measured the conditional, 2-point probability density 
P{yi\xi; r) from random realizations of the potential field. 
Results are presented in Fig. 2] (symbols) for a separation 
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in the range 2 < r < 2.5 h ^Mpc, and for two particu- 



lar values of 



± 2tT, where aj 



what follows. These measurements are compared to the 
following Gaussian conditional distribution : 



Piyt\xi;r) 



V 



(50) 



X exp 



0) 



2a? 



where, for shorthand convenience, we have introduced 
the variable Xi = ~ (xi). This conditional proba- 
bility is shown as the dashed curves. Figure [3] nicely 
demonstrates the validity of the Gaussian approxima- 
tion down to scales of the order of the grid resolution 
(about 1 /i^^Mpc). We also note that, in the distribu- 
tion P{yi\xi;r), the skewness increases by ~ 50 percent 
to reach « 0.089 at separation r ~ 2 /i~^Mpc whereas in 
P(y2\x2',r), the skewness is approximately 0.016 at the 
same distance. 

Thus far, we have restricted the comparison to the 
case in which the shear at Lagrangian positions and 
^2 is smoothed at the same length Rs- To further as- 
sess the validity of the Gaussian approximation in the 
limit d^ii « af, we take advantage of the fact that the 
joint distribution of the shear eigenvalues, evaluated at 
a single Lagrangian position Qi — snd two different 
smoothing lengths i?i and R2, also has the functional 
form of Eq. with = 7 [28]. Therefore, despite 

the limitation arising from the finite grid spacing, we 
can nevertheless probe the strongly correlated regime by 
studying the conditional probability density P(yi\xi;0) 
in the limit i?i « i?2- For illustration, let us consider 
the largest eigenvalue. We set i?i = 2 /i^^Mpc and take 
i?2 = 2.5, 3 and 5 h'^Mpc. The cross correlation at 
two different smoothing scales, computed following the 
procedure outlined in Sec. llVi yields af = 0.143, 0.136, 
and 0.108 (5l|, respectively. These values are used for the 
evaluation of the conditional density (jSTjl . As seen from 
Fig. [5l there is, as before, a remarkably good agreement 
between P(yi\xi;0) measured from random realizations 
of the potential $ and the conditional bivariate Gaus- 
sian equation (jSip , except for the very tails of the distri- 
butions. This strongly suggests that Gaussian statistics 
are an accurate approximation to the statistics of shear 
eigenvalues at all separations and smoothing scales. 
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FIG. 6: Same as Fig. [5] but showing the conditional proba- 
bility density P{diyi\dixi;0). Deviations from the bivariate 
Gaussian are more pronounced here than in the conditional 
distribution P{yi\xi;0). 



covariance tensor of the form 



(51) 



When Cii is known analytically, exact expressions can 
be derived fairly easily using the relation (diXidjyi) = 
—didjCii{r). For instance, the approximation p3p gives 



{diXidjyi) 



(52) 



where c = 87/270 — 9/107r and a prime denotes a deriva- 
tive with respect to r. Let us introduce the dimensionless 
variable diXi = (uq/ <Ti)diXi, where the g,- generally are 
the spectral moments of the density field [23| . 



(53) 



B. Gradients of the shear eigenvalues 



Using Eq. ([55]) as an analytic estimate of the correlation 
{diXidiyij , we find a variance 



The derivatives of any Gaussian random field X{q) 
with respect to the coordinate q also are Gaussian ran- 
dom fields since the differential operators d/dqi are lin- 
ear. To ascertain how much the spatial derivatives of the 
shear eigenvalues deviate from Gaussianity, we have ex- 
amined a number of correlation functions and conditional 
2-point statistics, focusing on the first derivatives diXi of 
the largest eigenvalue. This isotropic vector field has a 



1 /34 9 



{diil) = o 77 - ^ ~ 6.087 X 10^2 



3 V 45 57r 



(54) 



somewhat 10 per cent smaller than the value of « 6.60 x 
10~^ measured from the random realizations with Rs in 
the range 2 — 5 /i~^Mpc. The agreement can be improved 
by adding higher order terms in the truncated expansion 
equation 



11 



The conditional probability P{diyi\diXi;0) is shown 
in Fig. [HI for three different smoothing lengths i?i « R2 
(as in Fig[5|) . Cross correlations coefficients and variances 
are computed from the random realizations. Clearly, al- 
though deviations from Gaussianity are more pronounced 
than in the conditional 2-point statistics of the eigenval- 
ues discussed above, the agreement is still reasonable. 
However, we have found that it worsens significantly for 
the second derivative. In spite of this limitation, it would 
be valuable to assess whether the Gaussian approxima- 
tion provides a reliable description of the statistics of 
extrema of the shear eigenvalues, since the latter play a 
particular role in nonspherical collapse models @. 

VI. CONCLUSION 

We have explored the statistical correlation that arises 
in Gaussian initial conditions between the properties of 
the Hnear tidal shear didj^ at two distinct positions, 
thereby extending the work of [H [13, [H, [23, (H]. In 
Sec. Illli using asymptotic expansions, we derived exact 
closed form expressions for the joint distribution of shear 
components and shear eigenvalues as a function of the La- 
grangian separation. These results were applied to study 
the large-distance asymptotics of the correlation function 
of the shear eigenvalues and the shear principal axes. In 
Sec. IIV[ we presented interpolation formulae that accu- 
rately match the large- and small-scale behavior of the 
correlation of shear eigenvalues measured from random 
realizations of the gravitational potential. We also found 
that the alignment of the shear principal axes of thresh- 
olded regions is insensitive to the threshold height. How- 
ever, restricting the correlation to regions where all three 
eigenvalues are positive, introduces a dependence on the 
threshold density, which manifests itself as a strong sup- 
pression of the alignment for peak height less than v ^ 1. 
We emphasize that all these correlations are proportional 
to the density correlation on large scale. 

In Sec. |Vl we showed that the 2-point statistics of the 
shear eigenvalues closely follow the Gaussian statistics 
regardless of the separation and the smoothing length. 
Although we have not formally established that Gaussian 
multivariates comply with measurements of the ri-point 
distributions of shear eigenvalues, we speculate that the 
Gaussian approximation also holds for these multipoint 
distributions. Under this assumption, it should be fairly 
straightforward to apply the techniques and results ob- 
tained for Gaussian density fields to the shear eigenval- 
ues. 

Gaussian statistics provides also a reasonable descrip- 
tion of a number of conditional probability densities in- 
volving first derivatives of the shear eigenvalues. Note, 
however, that the agreement worsens noticeably for the 
second derivatives. This caveat notwithstanding, a Gaus- 
sian approximation should be adequate to understand, at 
least qualitatively, the clustering of extrema of the shear 
eigenvalues for instance. The mathematical framework 
laid down by [HI appears well suited for such a study. 



Our results can also be applied to the description 
of large-scale structures using the cosmic web approach 
based on the ellipsoidal collapse d, [13] ■ In light of our 
analysis, the conditional multivariate Gaussian describ- 
ing the joint distribution of the density, displacement 
field, and shear could easily be written down. As rec- 
ognized in Q, these statistics will prove useful for quan- 
tifying the properties of the mildly nonlinear fluctuations, 
which evolve into the network of clusters, filaments, and 
walls observed in the recent 2dF and SDSS galaxy surveys 
[l[ . Constraints on the tidal shear could also be included 
in topological measures such as Minkowski functionals, in 
an attempt to study the effect of nonspherical infall on 
the morphology of the primeval large-scale structures. 
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APPENDIX A: AVERAGING OVER THE 
RELATIVE ORIENTATIONS 

The integral (|22p over the relative orientation of the 
principal axis frames can be expressed as a hypergeo- 
metric function with the symmetric 3x3 matrices f3Ki 
and A2 as argument, 

p(3)^«A A ) - V V C.(Ai)C.(A2) 

(/3Ai,A2)-2^-2^ — . (Al) 

/C=0 K ^ ' 

Here designates summation over all partitions k h 
k of k, namely, over the ordered sequence of integer 
(fci, k2, ■ ■ ■ , kn) such that ki > k2 > ■ ■ ■ > kn > and 
^ fci = k. C„(X) are the zonal polynomials of the ma- 
trix X. They satisfy the relation (trX)*^ = EkC'«(X). 
We emphasize that, despite the use of matrix notation 
here, Ck(X) is a function of the eigenvalues of X solely 
and could thus be written as C'i^{x) for example. 

The zonal polynomial can be expressed in terms of the 
monomial symmetric functions 771^(2;). When k = 2 for 
instance, there are two zonal polynomials corresponding 
to the partitions (2) and (1,1) of 2, C(2) = m(2)(a;) -I- 
2/3to(i i)(a;) and C(i 1) = 4/3to(i j), where 

f \ 2,2,2 
m(2)(x) = Xi+X2+X^ 

W(i,i)(x) = a;ia;2 -I- si^s -f a;2a;3 . (A2) 

The value of these zonal polynomials at I is C(2) (I) — 5 
and C(i i)(I) = 4. There is a recurrence relation between 
the coefficients of 771^(2^) that determines C'k( X) uniquely 
once the coefficient of r77(fe) is given [s^, [11]. Note also 
that the functions mK{x) can be written in terms of the 
traces of power of X, trX*"' with k — 0, 1, ... . We refer 
the reader to [13] for further details. 
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